/*
** Modified by Neal May/June 2011
** Remodified by Neal Sept 2014 for Letters piece
** Stripped out lots of things to just compare Logit and clogit for standard normal x
** and fixed earlier version to compute RMS error and bias
** This program is a modified version of
** 
** Bias in Conditional and Unconditional Maximum Likelihood Estimation
**  Ethan Katz (ekatz@fas.harvard.edu)
**  last modified 02/06/01
** of Tom Coupe's
** date: 02 02 2005
** this modification solves the error in Katz's original program as it does allow to include fixed effects
** Tom Coup?  tcoupe@eerc.kiev.ua
** NB - 2014 The modification is now almost 100% so I am no longer calling this a modification 
** Modifed July 2015 to produce Tables 3 and 4
*/

version 14
clear all
/*This is done once and for all including setting any parms used and controlling looping*/
set matsize 2000
set seed 12947
set more off
capture log close
log using marginals.log, replace
/*the stuff below may change for different sets of runs*/
   matrix times =  (5\7\10\20\30\50\75\100\3)
local ntimes=rowsof(times)
matrix groups = (20\50\100)
local ngroups=rowsof(groups)
matrix conterm = (0)
matrix corrterm =(2)
local ncorr = rowsof(corrterm)
  local nnovar=0 /*number of units where all y are 0*/
local  sims=	1000 /* number of simulations */
 matrix xco=(1)
local nx = rowsof(xco)
local ncons=rowsof(conterm)
local rowres=`ncons'*`ntimes'*`ngroups'*`ncorr'*`nx'
matrix define resnew=J(`rowres',6,0)

matrix colnames resnew = N T RMSLOGIT   RMSCLOGIT RMSOLS  CORRTERM
/* 3 and 4 are RMS errors, 5 and 6 are signed errors, 7 is number of groups dropped */
noisily display "X is sum of SN mutliplied by xfecorr and the N fixed effects"
noisily display    " SIMS " `sims'
/*nothing changed below*/
local counter=0
  /* Now we loop over constants and T */

  forvalues xfn = 1/`nx' {
  forvalues consum=1/`ncorr' {
  forvalues groupnum = 1/`ngroups' {
forvalues connum=1/`ncons' {
forvalues  timenum = 1/`ntimes'  {
quietly {
local counter=`counter'+1
drop _all
set obs 10000
local xcoeff = xco[`xfn',1]
local xfecorr=corrterm[`consum',1]
local  time=times[`timenum',1]
 local n=groups[`groupnum',1]
local concoef= conterm[`connum',1] /* coefficient on constant term */
noisily display "sim number  " `counter' " n  " `n'  " t " `time'
local t=`time'			/* number of time periods */
local nomit=`nnovar'*`t'

local nt=`n'*`t'		
local n1=`n'+1
local n2=`n'+2
local t1=`t'+1


/*now we are going to generate the period and group variable*/
  
gen group = ceil(_n/`t') if _n<=`nt'
gen time = mod(_n-1,`t')+1 if _n<=`nt'
gen meanrsq = . /*rsq of x on fe1*/
gen errlogit = . /*absolute error in b for X coef of logit*/
    gen errlogittrue = . /*absolute error in b for X coef of logit w true dummy*/
gen errclogit = . /* absolute error in b for x coef of clogit*/
    gen errreg=.
    local i=1
while `i'<=`sims' {
/* now we are going to draw the true effects which do not vary by group */
generate fe=rnormal() if _n<=`nt'
generate fe1=fe
replace fe1 = fe1[_n-1] if group == group[_n-1]
/*we are going to give x some kurtosis by using a mild lognormal */
generate xstart= rnormal() if _n<=`nt'
 /*this is to get x corr with fe's* and then standardized */
egen x=std(`xfecorr'*xstart +fe1) if _n<=`nt'


/*this is to get the R^2 of x on fe1 */
  regress x fe1 in 1/`nt'
replace meanrsq=e(r2) in `i'
/*uncomment below if want to drop some groups with all y zero */
/*sort fefe group
local nomit=`n'*`nnovar'
replace fefe=-5 if _n<=`nomit' */

	/* draws Y from a Bernoulli distribution and estimates parameters */
gen p=	1/(1+exp(-(`concoef'+fe1+`xcoeff'*x))) if _n<= `nt'
                                           
                       
                     /*here is where we pick up the old sim code */
                       
	gen temp=uniform() if _n<=`nt'
	gen y=0
replace y=1 if temp<p



                    
capture clogit y x in 1/`nt', group(group) iterate(10)
           
                        if _rc!=0 | e(ic)>9 {
                            noisily  display "bad clogit"
                            
                            capture drop temp   y   p   x fe fe1 xstart
                         continue
                       }
      
local tc=_b[x]
         capture drop samp
         gen samp=e(sample)

capture 	logit y i.group  x if samp==1, iterate(10)
       
if _rc!=0 | e(ic)>9  {
                         noisily display "bad logit"
                         capture drop temp  y   p   x fe fe1 xstart samp
                         continue
                     }

/* predict pi and compute marginals over group and average*/

    capture drop predpi
predict predpi if samp==1
capture drop marg
gen marg=predpi*(1-predpi)*_b[x] if samp==1

capture drop truemarg
gen truemarg = p*(1-p)*`xcoeff' if samp==1

summar marg if samp==1
local predmarg=r(mean) 
summar truemarg if samp==1
local tmarg=r(mean)
       replace errlogit=(`tmarg'-`predmarg')^2 in `i'
       constraint 1 x=`tc' 
capture logit y i.group x if samp==1, constraint(1) iterate(10)
         if _rc!=0 | e(ic)>9  {
                         noisily display "bad logit"
                         capture drop temp  y   p   x fe fe1 xstart samp
                         continue
                     }
         capture drop predpicl
                     predict predpicl if samp==1
       capture drop margl
       gen margl=predpicl*(1-predpicl)*_b[x] if samp==1
       summar margl if samp==1
       local predmargl=r(mean)
replace errclogit=(`tmarg'-`predmargl')^2 in `i'

capture reg y x i.group if samp==1
    if _rc!=0   {
                         noisily display "bad logit"
                         capture drop temp  y   p   x fe fe1 xstart samp
                         continue
                     }
replace errreg=(_b[x]-`tmarg')^2 in `i'

capture drop temp   y   p   x fe fe1 xstart
	local i=`i'+1

      }
           
/*this is the end of a single sim */



/*done after all the sims for a given T and concoef are done */

    matrix resnew[`counter',1]=`n'  /*number of groups*/
        summar meanrsq
  matrix resnew[`counter',6]=r(mean) /*corrterm*/
  

matrix resnew[`counter',2]=`time' /*number of time periods*/

 





summarize errlogit
matrix resnew[`counter',3]=sqrt(r(mean))
summarize errclogit
matrix resnew[`counter',4]=sqrt(r(mean))
summarize errreg
matrix resnew[`counter', 5] = sqrt(r(mean))

noisily {
 
matrix list  resnew
}

}  /*quietly*/
}  /*xcoef*/
} /* corr*/
} /* groups */
} /* constant */
} /* time */
/*Done once at the end after all looping is done*/
noisily {
 
matrix list  resnew
}

capture erase logitrms2.csv
mat2csv, matrix(resnew) saving(logitrms2)

